package fem2.element;

import java.util.ArrayList;

import math2.MathUtilities;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.math.LinAlg;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;

/**
 * embedded crack element with constant opening crack
 * 
 * @author hbui
 * 
 */
public abstract class EmbeddedCrackElement extends StructuralElement {

	protected double[] uHatLocal; // buffer uHatLocal

	private boolean isActivated = true;

	private boolean isCrack;

	protected double crackLength;
	protected double area;
	protected double[] normal;
	protected double[] tangent;

	protected ArrayList<double[]> X = new ArrayList<double[]>();

	protected ArrayList<Integer> posDomain = new ArrayList<Integer>();
	protected ArrayList<Integer> negDomain = new ArrayList<Integer>();

	protected double ft; // tensile strength

	public EmbeddedCrackElement(MeshPart mp, MaterialModel mm, double ft) {
		super(mp, mm);
		this.ft = ft;
	}

	/**
	 * Compute enhanced strain operator using in compute sigma. sigma = Bu +
	 * G(xi) * alpha
	 * 
	 * @param xi
	 * @return the enhanced strain operator
	 */
	protected abstract IMatrix computeEnhancedStrainOperator(double... xi);

	/**
	 * Compute the projection operator for sigma in computing the local
	 * equilibrium. rlocal = int(G(x)^T * sigma, \Omega) + int(t, \Gamma)
	 * 
	 * @param xi
	 * @return the projection operator
	 */
	protected abstract IMatrix computeProjectionOperator(double... xi);

	/**
	 * get the parameters associated with the crack for the model
	 * 
	 * @return the crack opening parameter
	 */
	public abstract IVector getCrackOpeningParameter();

	/**
	 * get a representative stress in element to perform spectral decomposition
	 */
	protected abstract IVector getRepresentativeStress();

	/**
	 * 
	 * @return true if element is activated<br>
	 *         false if not
	 */
	public boolean isActivated() {
		return isActivated;
	}

	/**
	 * activate the element
	 * 
	 * @param active
	 */
	public void activate(boolean active) {
		this.isActivated = active;
	}

	/**
	 * 
	 * @return true if element is cracked<br>
	 *         false if not
	 */
	public boolean isCrack() {
		return isCrack;
	}

	/**
	 * 
	 * @return the normal vector pof crack
	 */
	public double[] getNormal() {
		return normal;
	}

	/**
	 * get list of nodes associated with the crack
	 * 
	 * @return
	 */
	public ArrayList<double[]> getX() {
		return X;
	}

	/**
	 * 
	 * @return the positive domain of the crack
	 */
	public ArrayList<Integer> getPositiveDomain() {
		return posDomain;
	}

	@Override
	public void Initialize() {
		super.Initialize();
		uHatLocal = new double[getMeshPart().countDofs()];
	}

	/**
	 * Typically, there are two possibilities for the normal direction. In
	 * Jirasek's triangle formulation, the normal always point to solitary node.
	 * However it is not necessary in Linder/Sancho method. Derive this method
	 * if you have other posibitites for normal direction <br>
	 * Note:<br>
	 * The normal direction will affect the positive domain
	 * 
	 * @param n
	 * @return the normal direction
	 */
	protected double[] selectNormal(double... n) {
		return n;
	}

	/**
	 * chec the onset criteria using Rankine criterion. Concurrently, compute
	 * normal, tangent, crack length,... when crack is initiated
	 */
	protected void checkCrackCriterion() {
		IVector sigma = getRepresentativeStress();
		double oxx = sigma.get(0);
		double oyy = sigma.get(1);
		double oxy = sigma.get(2);
		double t1 = (oxx + oyy) / 2;
		double t2 = Math.sqrt(Math.pow((oxx - oyy) / 2, 2) + Math.pow(oxy, 2));
		double o1 = t1 + t2;

		if (o1 >= ft) {
			isCrack = true;
			MeshPart mp = getMeshPart();

			/*
			 * compute crack center
			 */
			double[] crackCenter = mp.computeCenter();

			/*
			 * realize crack normal
			 */
			double n1, n2;
			if (Math.abs(oxy) < MathUtilities.tol14) {
				if (oxx > oyy) {
					n1 = 1;
					n2 = 0;
				} else {
					n1 = 0;
					n2 = 1;
				}
			} else {
				double g = Math.sqrt(Math.pow(oxx - o1, 2) + Math.pow(oxy, 2));
				n1 = -oxy / g;
				n2 = (oxx - o1) / g;
			}
			normal = selectNormal(n1, n2);
			tangent = new double[] { -normal[1], normal[0] };

			/*
			 * realize postive domain
			 */
			Node[] nodes = mp.getNodes();
			int nn = mp.countModes();
			// seggregate nodes on one side of the crack
			double a = normal[0];
			double b = normal[1];
			double c = -(a * crackCenter[0] + b * crackCenter[1]);
			for (int i = 0; i < nn; i++) {
				double[] p = nodes[i].getPosition();
				if (a * p[0] + b * p[1] + c > 0) {
					posDomain.add(i);
				} else {
					negDomain.add(i);
				}
			}
			// reverse if PosDomain is not in the direction of normal
			double[] p = nodes[posDomain.get(0)].getPosition();
			if (normal[0] * (p[0] - crackCenter[0]) + normal[1] * (p[1] - crackCenter[1]) < 0) {
				@SuppressWarnings("unchecked")
				ArrayList<Integer> tmp = (ArrayList<Integer>) posDomain.clone();
				posDomain = negDomain;
				negDomain = tmp;
			}

			/*
			 * compute crack length and two intersections
			 */
			// realize two edges intersect with crack
			int[] edge1, edge2;
			if (nn == 3) {
				if (posDomain.size() == 1) {
					edge1 = new int[] { posDomain.get(0), negDomain.get(0) };
					edge2 = new int[] { posDomain.get(0), negDomain.get(1) };
				} else if (posDomain.size() == 2) {
					edge1 = new int[] { negDomain.get(0), posDomain.get(0) };
					edge2 = new int[] { negDomain.get(0), posDomain.get(1) };
				} else {
					throw new Error("Something wrong with the domain");
				}
			} else if (nn == 4) {
				if (posDomain.size() == 1) {
					int cnt = 0;
					int[] adjNodes = new int[2];
					for (int i = 0; i < 3; i++) {
						if (Math.abs(posDomain.get(0) - negDomain.get(i)) != 2) {
							adjNodes[cnt++] = negDomain.get(i);
						}
					}
					edge1 = new int[] { posDomain.get(0), adjNodes[0] };
					edge2 = new int[] { posDomain.get(0), adjNodes[1] };
				} else if (posDomain.size() == 2) {
					if (Math.abs(posDomain.get(0) - negDomain.get(0)) != 2) {
						edge1 = new int[] { posDomain.get(0), negDomain.get(0) };
						edge2 = new int[] { posDomain.get(1), negDomain.get(1) };
					} else {
						edge1 = new int[] { posDomain.get(0), negDomain.get(1) };
						edge2 = new int[] { posDomain.get(1), negDomain.get(0) };
					}
				} else if (posDomain.size() == 3) {
					int cnt = 0;
					int[] adjNodes = new int[2];
					for (int i = 0; i < 3; i++) {
						if (Math.abs(negDomain.get(0) - posDomain.get(i)) != 2) {
							adjNodes[cnt++] = posDomain.get(i);
						}
					}
					edge1 = new int[] { negDomain.get(0), adjNodes[0] };
					edge2 = new int[] { negDomain.get(0), adjNodes[1] };
				} else {
					throw new Error("Something wrong with the domain");
				}
			} else {
				throw new Error("invalid number of nodes");
			}
			// compute intersection with edge 1
			double[] N1 = nodes[edge1[0]].getPosition();
			double[] N2 = nodes[edge1[1]].getPosition();
			double xi1 = N1[0], yi1 = N1[1];
			double xi2 = N2[0], yi2 = N2[1];
			double a1 = yi2 - yi1;
			double b1 = xi1 - xi2;
			double c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] X1 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			X.add(X1);
			// compute intersection with edge 2
			N1 = nodes[edge2[0]].getPosition();
			N2 = nodes[edge2[1]].getPosition();
			xi1 = N1[0];
			yi1 = N1[1];
			xi2 = N2[0];
			yi2 = N2[1];
			a1 = yi2 - yi1;
			b1 = xi1 - xi2;
			c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] X2 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			X.add(X2);
			// compute crack length
			crackLength = LinAlg.diffNorm2(2, X1, X2);

			/*
			 * compute area
			 */
			area = getFace().computeArea();
		}
	}

}
